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Convective  flow  with  subcritical  instability 

John  C.  Neu 

Department  of  Mathematics.  Stanford  University,  Stanford.  California  94305 
(Received  4  April  1980;  accepted  - '  October  19811 

An  asymptotic  analysis  of  subcritical  instability  in  double  diffusive  convection  is  presented.  Using  a  modified 
perturbation  method,  a  Landau  equation  that  determines  how  the  amplitude  of  the  convection  evolves  in  time 
is  derived  From  the  Landau  equation,  it  is  found  that  in  certain  cases,  stable  finite  amplitude  convection  can 
exist  even  when  the  rest  state  with  no  flow  is  locally  stable.  The  perturbation  analysis  complements  and 
unifies  previous  work  which  is  primarily  qualitative  or  numerical  in  character. 


I.  INTRODUCTION 

There  are  fluid  flows  which  are  stable  against  infinite¬ 
simal  perturbations,  yet  finite  amplitude  disturbances, 
once  created,  may  persist  and  even  grow.  We  state 
that  such  flows  exhibit  subcritical  instability.  Certain 
geophysical  and  astrophysical  convection  processes 
provide  prime  examples.  Veronis1  discovered  subcri¬ 
tical  instability  in  the  process  of  double  diffusive  con¬ 
vection.  Here,  the  basic  thermal  instability  is  mod¬ 
ified  by  the  presence  of  a  solute  which  introduces  a 
second  bouyancy  force  in  addition  to  the  usual  effect  of 
thermal  expansion.  Previous  researchers  into  this 
subject  have  proceeded  along  the  lines  of  a  simplified 
semi-quantitative  description,  or  a  detailed  numerical 
simulation.  In  this  paper,  we  present  a  quantitative 
asymptotic  analysis  of  double  diffusive  convection.  Cer¬ 
tain  results  of  previous  analyses  are  consequences  of 
the  unified  treatment  presented  here. 

To  pose  the  proper  goals  of  the  analysis,  we  consider 
the  physical  mechanism  of  subcritical  instability  in 
double  diffusive  convection.  In  a  fluid  with  no  solute, 
where  purely  thermal  convection  takes  place,  the  ex¬ 
istence  of  convection  due  to  an  adverse  tempera¬ 
ture  gradient  is  determined  by  a  single  dimen¬ 
sionless  parameter  R  called  the  Rayleigh  number.  R 
is  proportional  to  A7'/v,  where  AT  is  the  temperature 
difference  across  the  fluid  layer  and  v  is  the  viscosity. 
There  is  a  critical  Rayleigh  number  R=R* ,  such  that 
the  static  solution  corresponding  to  heat  transfer  purely 
by  thermal  convection  without  fluid  flow  is  stable  for 
R<R*  and  unstable  for  R>  R*.  For  R  >  R* ,  the  stable  state 
is  steady  convection.  The  amplitude  A  of  the  convection,  as 
measured  by  the  velocity  at  a  fixed  point  in  the  flow, 
increases  like(R  -  fl*)‘/2for  0<  ft  -  R*«  1.  Curve  (i)  in 
Fig.  1  is  a  plot  of  the  amplitude  as  a  function  of  Rayleigh 
number  for  purely  thermal  convection  and  is  typical 
of  systems  that  exhibit  supercritical  bifurcation. 

We  now  consider  the  modifications  due  to  the  presence 
of  a  solute.  We  assume  that  a  higher  solute  concentra¬ 
tion  is  maintained  on  the  bottom  of  the  fluid  layer  than 
on  the  top,  so  that  the  solute  appears  as  a  stabilizing 
agent.  Hence,  the  critical  Rayleigh  number  increases  to 
a  value  ft,^ft*.  We  furthei*  assume  that  the  solute 
diffusivity  is  sufficiently  small  so  that  the  mixing  due  to 
convective  motions  is  not  overcome  by  the  effect  of 
fixed  boundary  concentrations.  Under  these  conditions, 
it  may  be  possible  for  finite  amplitude  convection 
to  exist  at  thermal  Rayleigh  numbers  ft<ft,, 
because  the  convective  motion  maintains  a  nearly  uni¬ 


form  solute  concentration  in  the  interior  of  the  roll,  and 
the  stabilizing  solute  gradient  does  not  appear  there. 
Curve  (ii)  in  Fig.  1  shows  the  conjectured  bifurcation 
diagram  for  double  diffusive  convection.  The  locus  of 
stable  steady  states  is  indicated  by  the  solid  lines,  and 
the  locus  of  unstable  steady  states  by  hatched  lines. 
According  to  this  picture,  there  is  stable,  finite  ampli¬ 
tude  convection  for  subcritical  Rayleigh  numbers  in  an 
interval  Rm<  R<  ft,.  We  state  that  the  convective  pro¬ 
cess  exhibits  subcritical  bifurcation.  The  goal  of  the 
present  work  is  to  provide  an  asymptotic  description  of 
this  phenomenon. 

Huppert  and  Moore2  discovered  a  stable  branch  of  sol¬ 
utions  that  represent  finite  amplitude  convection.  These 
solutions  were  computed  both  numerically  and  by  analy¬ 
tic  approximation  based  on  mean  field  equations.  This 
branch  is  represented  by  the  solid  portion  of  curve  (ii) 
in  Fig.  1.  In  addition,  they  discovered  a  branch  of  un¬ 
stable  solutions  emanating  from  a  bifurcation  point  at 
R-Rs,S  =  0.  This  branch  is  represented  by  the  hatched 
portion  of  curve  (ii).  In  this  paper,  we  verify  the  con¬ 
jecture  of  Huppert  and  Moore  that  the  stable  and  un¬ 
stable  branches  of  solutions  connect  each  other  in  a 
continuous  fashion  to  form  a  single  branch  of  solutions. 
In  the  limiting  case  with  small  solute  diffusivity  and 
small  solute  gradients,  we  present  a  unified  analysis 
which  describes  the  stable  and  unstable  branches,  and 
their  joining  to  form  one  continuous  family  of  solutions. 

In  Sec.  II,  we  give  the  mathematical  formulation  of  (he 
double  diffusive  convection  problem.  In  Sec.  Ill,  we 
present  the  results  of  linearized  stability  theory  for 
periodic  rolls.  In  Sec.  IV,  we  obtain  a  first  qualitative 
description  of  the  subcritical  instability  by  means  of  an 
averaging  procedure  due  to  Stuart.  His  procedure,  al¬ 
though  qualitative,  nevertheless  gives  valuable  insight 
into  the  conditions  under  which  subcritical  instability 
occurs.  In  particular ,  it  suggests  the  proper  limiting 
case  for  an  asymptotic  analysis.  In  Sec.  V,  we  begin 
a  formal  asymptotic  analysis  analogous  to  the  work  of 
Keller  and  Kogelman  on  the  Benard  problem.2  The 
eventual  result  of  this  analysis  is  a  Landau  equation  that 
describes  the  slow  temporal  variation  of  spatially  per¬ 
iodic  rolls.  The  derivation  of  the  Landau  equation  in¬ 
volves  an  elliptic  boundary  value  problem  for  the  solute 
concentration.  The  completion  of  the  analysis  does  not 
require  the  solution  of  the  elliptic  problem,  but  rather 
a  certain  inner  product  of  the  solution.  There  is  a  very 
convenient  variational  formulation  suggested  by  Keller'' 
for  computing  such  inner  products.  In  Sec.  VI,  we 
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FIG.  1.  Amplitude  of  convection  A  as  a  function  of  Rayleigh 
number  R. 


The  condition  w  =  0  at  z  =  0  and  *  =  1  expresses  the  con¬ 
finement  of  the  fluid  between  the  planes  2  =  0  and  2  =  1. 
u,=  0  is  the  “free  surface”  condition  of  no  tangential 
stress.  The  conditions  (3b)  follow  from  the  choice  of 
fixing  the  temperature  and  solute  concentration  at  the 
boundaries. 

The  system  (2a)-(2c),  (3a),  (3b)  is  the  mathematical 
formulation  of  a  double  diffusive  convective  process. 
The  problem  is  to  solve  Eqs.  (2a)-(2c)  for  9,  and  £ 
subject  to  the  boundary  conditions  (3a)  and  (3b). 


present  the  variational  procedure  and  its  application  to 
our  analysis.  Finally,  in  Sec.  VII  we  present  the  com¬ 
pleted  Landau  equation  and  its  predictions  concern  sub- 
critical  instability. 

II.  THE  EQUATIONS  OF  MOTION 

We  consider  a  convective  process  taking  place  in  a 
fluid  layer  of  thickness  d  with  both  surfaces  free.  Gra¬ 
vity  acts  uniformly  with  magnitude  g  in  the  z  direction. 
The  density  of  the  fluid  depends  on  the  temperature  T 
and  solute  concentration  S  according  to  an  equation  of 
state 

p  =  Po[l-a,(T-ro)  +  a#(S-S0)].  (1) 

Here,  p0  is  the  density  when  T  =  T0  and  S=S0;  a,  and 
a3  are  both  positive  constants.  Constant  temperatures 
and  solute  concentrations  are  maintained  on  the  free 
surfaces,  with  the  higher  values  on  the  bottom  surface. 
We  denote  the  differences  in  temperature  and  solute 
concentration  from  bottom  to  top  by  AT  and  AS.  In  this 
analysis,  we  consider  two-dimensional  motions  which 
depend  only  on  x  and  z.  We  nondimensionalize  by  mea¬ 
suring  all  lengths  in  units  of  the  fluid  layer  thickness 
d,  and  time  in  units  of  d2/nt,  where  k,  is  the  thermal 
diffusivity.  Under  these  conditions,  the  dimensionless 
equations  of  the  Boussinesq  approximations  are 

(8,  -  ffVl)VJi p  =  -oR6t  +  oSZt  -  ip,V2ip,  +  <p,V2<p.  ,  (2a) 
U,  -  v2)e  +  if,,  =  -4>tet  +  ,  (2b) 

(d'-TV2)  £  +  *,=  -*,£„ +«-£,,  (2c) 

for  0<  2<  1.  Here,  ip  is  the  stream  function  from  which 
one  may  determine  the  horizontal  and  vertical  velo¬ 
cities  according  to  u  =  <pt  and  w  =  -ipx.  6  and  2  are  the 
deviations  of  temperature  and  solute  concentration  from 
the  rest  state  with  no  fluid  flow,  ip,  9,  and  £  are  the 
unknowns  to  be  solved  for.  There  are  four  dimension¬ 
less  numbers:  o,  the  Prandtl  number,  is  the  ratio  of 
viscous  diffusivity  v  to  thermal  diffusivity  k,.  t  is  the 
ratio  of  solute  diffusivity  k,  to  thermal  diffusivity.  R 
sgat^Td2/K,v  and  S  *gat&Sif/ictv  are  the  thermal  and 
solute  Rayleigh  numbers,  respectively.  Equation  (2a) 
governs  the  momentum.  The  bouyancy  force  terms  are 
-oR9x  and  <rSZx.  Equations  (2b)  and  (2c)  express  the 
conduction  of  heat  and  solute. 

The  boundary  conditions  at  the  free  surfaces  are 
ip  =  -  0 ,  on  2  =  0,1.  (3a) 

|  0  =  £  =0,  on  2  =  0,1  .  (3b) 

The  conditions  (3a)  imply  w  =  ~ipx  =  0  and  u,  =  ipu  =  0. 


III.  LINEAR  STABILITY  THEORY 

The  rest  state  with  no  convective  motion  is  charac¬ 
terized  by  ip,  9 ,  £  =0.  From  a  linearized  analysis  of 
the  problem  (2a)-(2c),(3a),(3b),  we  see  that  this  rest 
state  is  unstable  against  spatially  periodic  perturba¬ 
tions  of  the  wavenumber  k  if  the  thermal  Rayleigh  num¬ 
ber  R  exceeds  the  critical  value 

RJ  =  R*  +  S/ t,  R*=(i2  +  k‘) 3/k2.  (4) 

R*  is  the  critical  value  for  purely  thermal  convection 
withS  =  0.  If  S>  0,  then  Rt>  R".  We  see  that  the  pre¬ 
sence  of  a  solute  whose  concentration  increases  with 
depth  has  a  stabilizing  effect.  If  R=RS,  we  have  the 
condition  of  neutral  stability  for  which  the  linearized 
problem  admits  a  time  independent  solution.  This  sol¬ 
ution  is  given  by 

ip  —  at  sinfcxsinirz,  9  =aecos*xsinjrz  , 

£  =  0ccosfeA'Sinsz ,  (5) 

where 

Of  ~-(k/n2  +  k2)at,  aL  =  -(l/r)(k/n2  +  k2)at  .  (6) 

Figure  2  depicts  the  streamlines,  which  are  curves  of 
constant  ip. 

IV.  THE  QUALITATIVE  THEORY  OF  SUBCRITICAL 
INSTABILITY 


Chandrasekhar5  has  analyzed  finite  amplitude  thermal 
convection  by  an  averaging  procedure  called  Stuart’s 
method.  We  generalize  his  work  to  the  case  of  double 
diffusive  convection.  The  essence  of  the  procedure  is 
simple.  From  the  equations  of  motion  (2a.) — (2c ),  (3a), 
(3b),  we  derive  three  integral  relations  for  steady  con¬ 
vection  which  express  the  balances  between  the  creation 
and  removal  of  mechanical  energy,  heat  and  solute. 
These  integral  relations  are  given  by 


/<*•  V2ip  }dz-R  j  (ipx9  )dz  +  S  j  {pxs)dz-  0  ,  (7a) 

f  <9V23>dz  -  f  (pt9)dz  =  f  (ipt9?dz  -if  {ipx9)dz\  , 

■'o  •'n  */o  \r  o  / 


(7b) 


FIG.  2.  Streamlines  for  neutrally  stable  convection. 
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rf  (sV2s)dz  -  f  (ip  s)dz 

=  ~  f  (<pts)~dz-^  j  {ipts)dzj  .  (7c) 

Here,  the  brackets  (  )  denote  the  average  over  a  hori¬ 
zontal  plane,  while  Safi-®)  and  saS-  <Z>. 

Near  marginal  stability  where  0 <R-Ra«  1,  the  con¬ 
vection  is  of  small  amplitude  and  the  field  variables  tp, 

8,  and  Z  are  well-approximated  by  the  values  (5)  of  the 
linearized  theory.  Upon  substituting  these  values  into 
the  integral  relations,  we  obtain  the  following  system  of 
equations  for  the  amplitudes  a0,  <%,  and  ac: 


Ur2  +  i2)a2  +  Rfeq,  at  -  Skavat  =  0  , 

(8a.) 

(k2  +  z2)<4  +  ka9at  =  -  a  k raja2  , 

(8b) 

r(k2  +  n2)a2c  +  kacat  =  -  (1/8 T)k2a2za\  . 

(8c) 

One  solution  of  these  equations  is  at  =  a9  =ac=0,  cor¬ 
responding  to  the  rest  state  with  no  convection;  but 
there  are  nonzero  solutions  as  well.  By  eliminating  a „ 
and  ac,  we  obtain  the  following  equation  relating  the 
streamfunction  amplitude  at  to  the  Rayleigh  numbers 
R  and  S: 

8(»2  +  *2)rj"2. 


(9) 

Figure  3  is  a  graph  based  on  Eq.  (9)  of  the  amplitude  a 
as  a  function  of  (R  -R’J/t2  for  various  values  of  S  and  r. 

If  S  a  0  and 

A3(1/t3-1/t)S-R*<0,  (10) 

then  (9)  yields  supercritical  bifurcation  represented  by 
curve  (i).  If  A>0,  there  is  subcritical  bifurcation  as 
represented  by  curve  (iii).  Curve  (ii)  represents  the 
intermediate  case  A  =  0,  in  which  the  bifurcation  curve 
makes  fourth-order  contact  with  its  vertical  tangent  at 
the  point  P. 

V.  PERTURBATION  ANALYSIS 

We  present  a  quantitative,  asymptotic  analysis  of 
double  diffusive  convection  to  compliment  previous 
works  in  the  subject  which  have  been  either  qualitative 
or  numerical  in  character.  The  asymptotic  limit  to  be 
considered  is  S -  0.  We  take  S  =  t  ,0<(  «  1.  From  the 
qualitative  theory  in  Sec.  IV,  we  expect  subcritical  bi¬ 
furcation  to  occur  when  A^d/r*  -  1/t)S_R*  >0,  Since 
S  is  of  order  « ,  we  see  that  subcritical  bifurcation  oc- 


F1G.  3.  Bifurcation  diagrams  based  on  Stuart's  method. 


curs  if  t,  the  ratio  of  solute  to  thermal  diffusivities,  is 
of  order  e'/3.  We  set  r  =  et/3T0,  where  tu  is  of  order 
unity.  Given  that  S  and  t  have  magnitudes  i  and  f'r\ 
the  proper  scaling  of  the  remaining  parameters  and 
variables  is  determined  from  the  following  considera¬ 
tions. 

(i)  The  fluxes  of  solute  due  to  diffusive  and  convec¬ 

tive  effects  should  balance.  In  the  solute  conduction 
equation  (2c)  the  term  tv2£  represents  the  effect  of 
diffusion,  and  the  terms  +  4  Z  represent  the  ef¬ 

fect  of  convection.  Since  t  is  0(r f/3),  we  see  that  tp  is 
0(f1/3). 

(ii)  For  purely  thermal  convection  with  0<R  -R* 

«  1,  the  amplitudes  of  ip  and  8  are  both  proportional  to 
(R  -R*)1/2.  Although  this  relation  is  not  valid  quantita¬ 
tively  for  double  diffusive  convection,  we  assume  that  it 
indicates  correct  orders  of  magnitude.  Hence,  8  is  also 
0(€1/3)  and  R -R*  is  0(t2/3).  For  Rayleigh  numbers 
R  with  R  —  R*  =  0(t 2/3),  the  growth  rate  of  temporal  in¬ 
stability  is  0(c2/3).  Hence,  the  time  scale  for  the 
growth  of  convection  is  0(c"2/3). 

(iii)  The  magnitude  of  the  thermal  bouyancy  force  due 
to  thermal  expansion  is  Re.  If  R  >  R* ,  a  solute  free 
fluid  experiences  convective  instability.  Hence,  we 
may  think  of  R*8  as  the  portion  of  the  thermal  bouyancy 
force  that  is  required  to  overcome  the  stabilizing  effect 
of  viscosity,  while  ( R~R*)9  represents  the  actual  mo¬ 
tive  force  of  the  instability.  In  the  case  we  consider, 
where  a  solute  is  present,  we  expect  that  this  destabil¬ 
izing  force  will  be  balanced  in  magnitude  by  the  solute 
bouyancy  force  SZ.  Since  ( R-R*)8  is  0[(e2/’)U 1  /3)] 

=  0{e  )  and  S  =  f ,  we  find  that  Z  is  0(1). 

On  the  basis  of  the  discussion  in  (i),  (ii),  and  (iii), 
we  adopt  the  following  scaling  of  the  variables: 

<p  =  e'n*,  8  =e  1/36,  Z=M,  (11) 

R  =  R*  +  e2/3r,  t=e2nT. 

We  seek  asymptotic  solutions  for*,  0,  and  ^  in  the  form 

e-ep  +  e1/3e‘  +  £2/10!+ . .  .  ,  02) 

*~S0+E1/3S‘  +  e2/3S2+  .  .  .  . 

We  discuss  the  leading  order  solution.  From  Eqs.  (2a) 
and  (2b)  we  find  that  ♦°  and  O'  satisfy 

V24'°  -R*0^  =  0,  V2e°-*°=0.  (13) 

The  solutions  which  satisfy  the  boundary  conditions  +° 

=  =  0,  9°  =  0  at  z  =  0  and  z  =  1  are 

#°  =  a  sin/rxsinnz,  0°=  -[k'(k2  +  i2)]a'cos/fexsin»z  . 

(14) 

Here,  a  =  a(T)  represents  the  time  varying  amplitude 
of  convection.  The  goal  of  this  analysis  is  to  find  its 
governing  equation. 

From  the  solute  conduction  Eq.  (2c)  we  find  that  £°, 
satisfies 


-In.  ^£0  _  ^"^0  »_  ^rv 


(15) 
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If  we  set  0o  =  a  sintfexsinjz,  this  becomes 

(T0/a)V25°-  (*>,£"+  <^,4^  =  <*>,.  0  =  sinfe.vsin!7z  .  (16) 

Due  to  the  periodicity  of  0  in  the  x  direction,  it  is  suf¬ 
ficient  to  solve  for  £°  inside  the  single  convection  cell 
depicted  in  Fig.  3  with  0  «  z  «  1  and  0  «  x  «  v/k.  The 
boundary  condition  on  2  =  0  and  2=1  is  5“  =  0.  Along  the 
interfaces  of  the  convection  cell  with  its  neighbors, 
symmetry  dictates  £“  =  0. 

To  complete  the  leading  order  description,  we  need 
the  evolution  equation  of  the  amplitude  at(T).  This  is 
found  in  the  process  of  solving  the  higher  order  equa¬ 
tions.  We  solve  the  equations  for  the  first  order  cor¬ 
rections  #l  and  0l  and  substitute  the  results  into  the 
equations  for  the  second  order  corrections  02  and  CF. 
The  second  order  equations  have  a  solvability  condition 
which  yields  the  governing  Landau  equation  for  a(T). 

The  analysis  is  straightforward  and  gives  the  result 

-~<fe2  +  *2)ar  =  +  ^  +  U\  *„> , 

Q  R  +  It  o 


where  (£n,  <pt )  is 

(17) 

4i>  /*l  /•'/* 

u°, 4>j=~rJ  f  e<t>tdxdz. 

*  •'o  */0 

(18) 

The  first  two  terms  on  the  right-hand  side  of  (17)  ap¬ 
pear  in  the  usual  Landau  equation  for  purely  thermal 
convection.  The  term  (£°,  0^)  represents  the  effect  of 
solute.  It  is  a  definite  function  of  the  amplitude  a  be¬ 
cause  a  appears  as  a  parameter  in  Eq.  (15).  A  varia¬ 
tional  principle  suggested  by  Keller4  provides  a  very 
convenient  tool  for  estimating  the  functional  dependence 
of  (S°,  0„)  on  a. 

VI.  A  VARIATIONAL  PRINCIPLE  AND  ITS 
APPLICATION 

Let  u,a,b  be  functions  defined  in  a  bounded  domain. 

We  wish  to  compute  the  value  of  the  inner  product 
(u,6)  when  u  is  the  solution  of  the  inhomogeneous  ellip¬ 
tic  equation 

Lu  =a  ,  (19) 

subject  to  homogeneous  boundary  conditions.  This 
problem  has  a  simple  variational  formulation:  Let  u* 
be  the  solution  of  (19)  and  let  v*  be  the  solution  of  the 
adjoint  problem 

L'v  =  b.  (20) 

We  define  a  functional  g(u,v)  by 

g(u,v)  =  (a,v)  +  (u,b)  ~{Lu,v) .  (21) 

A  simple  calculation  shows  that 

g(u*  +  oltv*  +  0j)  =  («*,6)-(L<71,<t2)  .  (22) 

We  see  that  g  attains  its  stationary  value  when  u  =n* 
and  v-v*,  and  that  the  stationary  value  is  precisely 
the  required  inner  product  (a*, 6).  In  the  self-adjoint 
case,  the  stationary  value  is  also  extremal.  Since  the 
problem  we  consider  is  not  self-adjoint,  the  stationary 
value  is  not  necessarily  extremal. 


The  application  to  the  analysis  of  Sec.  V  is  clear. 

We  apply  variational  principle  with 

L  =  (T0/a)V2  -<t>tdt  +  <ptdt,  a  =  b  =  <px .  (23) 

The  functional  g(u,v)  is 

g(u,v)  =  (<px,  v)  +  (a,  0,)  -  [<T/2)V2a  -  <Pjit  +  , v ]  , 

where  the  inner  product  is  the  one  defined  by  (18). 

To  compute  the  stationary  value  of  g,  it  appears  that 
one  must  perform  variations  with  respect  to  both  a  and 
ti;  but,  there  is  a  symmetry  in  our  particular  problem 
which  allows  a  simplification.  The  adjoint  operator  is 
V  =  (  tJ  a )  V2  +  9,0,  -  a,  0,  =  (  tJ  a )  V2  +  0,  a,  -  <bx  a, . 

,  (25) 

Hence,  the  adjoint  equation  Uv  —  b  is 

(t J  a)V2v  +  <ptvt  -  <t>xv,  =  <px  .  (26) 

The  function  0(x,z)  =  sinfeAsinirz  has  the  property 
0(x,z)  =  0(x,  1  -z).  As  a  consequence  of  the  symmetry, 
we  find  that  the  adjoint  Eq.  (26)  becomes  equivalent  to 
the  original  equation  Lu  =a  under  the  change  of  variable 
2-  2  =  1-  2.  Hence,  the  solution  v*  of  the  adjoint  prob¬ 
lem  is  obtained  by  substituting  1-2  for  z  in  a*(*,z); 
that  is, 

v*(x,  z)  =  u*(x,  1  -  z) .  (27) 

As  a  result,  the  stationary  value  of  g(u,v)  is  equal  to  the 
stationary  value  of 

g(u)=g[u(x,z),u(x,l  -z)].  (28) 

To  estimate  the  stationary  value  of  g(u),  we  substitute 
into  (27)  a  simple  approximation  to  u  with  undetermined 
parameters  and  then  compute  the  stationary  value  with 
respect  to  the  parameters.  To  third  order  in  ar/r0,  the 
solution  of 

Lu  =  (Tq/ad^u  -  <t>,ux  +  =  0,  =  a  ,  (29) 

subject  to  boundary  conditions 

«  =  0  on  j  =  0,1,  ux  =  0  on  x=0,f/k  ,  (30) 

has  the  form 

u  —  Ci  cos  kx  sin  iz  +  c2  sin2iz  +  c3  coskx  sin  3zz  . 

(31) 

This  is  a  natural  choice  for  u.  We  substitute  (31)  into 
the  expression  (28)  for  g.  The  stationary  value  of  g 
with  respect  to  the  parameters  c,,  c2,  and  c3  is 


**  =  -v$*(*’  +  ,2)>/2A  1+3A2  +  2^^A4 


x  (1+2A2)'2  , 

(32) 

where 

,  k  a 

A  *  [3(fc2  +  f2)]1  12  T„  ' 

(33) 

We  take  g*  as  our  approximation  to  the  inner  product 

The  accuracy  of  the  approximation  is  easily  accessed. 
The  variational  procedure  generates  values  of  c,,c2, 
c,,  so  that  the  resulting  expression  (31)  for  u  is  correct 
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(34) 


to  third  order  in  A.  That  is, 

/<=«*  +  < 7, 

where  u *  is  the  exact  solution,  and  a  is  0(A4). 
(22)  we  find  for  this  value  of  u  that 


From 


g*  =g(u)  =  g[u(x,z),u(.X,  1  -  2>] 

=  (j<*,b)-[<7(*,z),<Hx,l  -*)].  (35) 

Since  the  operator  L  in  (29)  is  0(\/A )  and  a  is  0(A4),  we 
find  that  the  error  term  [Lo(x,z),<7U,l  -z)]  is  O (A7). 
Hence,  the  approximation  of  (2)°,  4>t)  in  (32)  has  error 
OW7). 


VII.  THE  LANDAU  EQUATION  AND  ITS 
CONSEQUENCES 


Using  the  estimate  of  determined  in  Sec.  VI, 

the  Landau  Eq.  (17)  can  be  written  as 


lAT=r'A-A3-nA 


l  +  3A2  +  2 


(k2  +  5  tr2\ 
\k2  +  9ir2)' 


x  (1  t  2A2y2  +OW), 


(36) 


where  x  =  (l  +  a)/[OT^(fez  +  j2)],  r'ar/U?*^)  and  4 
sl/fl’rj.  Steady  rolls  have  A r  =  0.  Their  amplitudes 
satisfy 


r'  =  A-+  n 


1  +  3A2  +  2 


\k2  +  9v2)  . 


x  (1+2 A2)’2  +  OiA*). 


(37) 


Figure  4  depicts  A  as  a  function  of  r‘  for  various  values 
of  4  and  k  =  t  The  curves  with  4  >  1  correspond  to 
suberitical  bifurcation.  On  any  given  curve,  the  por¬ 
tions  with  dr'  dA2  >  0  correspond  to  stable  steady 
states,  while  the  portions  with  dr'‘dA2<  0  correspond 
to  unstable  steady  states.  The  unstable  portions  are 
indicated  in  Fig.  4  by  hatched  curves. 

We  consider  the  asymptotic  validity  of  the  bifurcation 
diagrams.  We  expand  (37)  as 


FIG.  5.  Comparison  between  Stuart’s  method  and  variational 
procedure. 

r'  =  4  -  (4  -  1W2  +  2 ti[(k2  +  5ir2)/(k2  +  9ir2)}44  +  OM6) . 

(38) 

The  characteristic  shape  of  the  suberitical  bifurcation 
curves  in  Fig.  4  is  due  to  the  balance  of  the  quadratic 
and  quartic  terms  in  (38).  Hence,  A2  is  0(4  -  l)2. 

The  error  V(. A1')  will  be  negligible  to  leading  order  if 
(4  -  1)4»  (4  -  l)6,  or  4  -  1«  1. 

We  compare  the  result  (37)  of  the  variational  analy¬ 
sis  with  the  result  (9)  obtained  by  Stuart’s  method. 
Setting  R  =R*  +€2/V  =  R*(l  +«2/1^r'),  a  =  A,  S  =  <,  and 
t  ~e'/3T0,  we  find  that  (9)  becomes 

r'=A2  +  4/(1  +  A2)  +  0(e2/3) .  (39) 

Figure  5  shows  graphs  of  A  vs  r’  for  4  =  2  as  obtained 
by  Stuart’s  method  and  the  variational  procedure. 


A 


FIG.  4.  Bifurcation  diagrams  based  on  variational  method. 


VIII.  COMPARISON  WITH  RELATED  WORK 

Most  numerical  solutions  of  the  double  diffusive  con¬ 
vection  problem  are  performed  at  values  of  the  thermal 
and  saline  Rayleigh  numbers  for  which  the  flow  is  fully 
nonlinear.  In  these  cases,  the  perturbation  analysis 
does  not  apply.  In  the  number ical  work  of  Huppert  and 
Moore,2  typical  values  of  4  =S/R’,r>  range  between  10 
and  1,000,  whereas  the  results  of  the  perturbation  analy¬ 
sis  presented  in  Sec.  VII  are  valid  for  l4-ll«l.  It 
seems  that  the  perturbation  method,  which  gives  accur¬ 
ate  analysis  locally,  does  not  provide  global  results. 

In  this  respect,  certain  semi-quantitative  methods  per¬ 
form  better  by  providing  crude,  but  globally  applicable 
results. 

Veronis1  provided  the  first  semi-quantitative  analysis- 
of  suberitical  convection  through  the  use  of  a  truncated 
modal  representation.  The  result  of  his  analysis  is  a 
fifth-order  system  of  ordinary  differential  equations  for 
the  amplitudes  of  the  fundamental  modes  and  second- 
harmonic  corrections.  From  these  equations,  he  de- 
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TABLE  1.  Critical  thermal  Rayleigh  number  for  finite  ampli¬ 
tude  convention  as  a  function  of  saline  Rayleigh  number. 


s 

Numerical 

^mln 

Formula  (41) 

Formula  (42) 

10s 

20.4 

18.6 

16.9 

103'1 

52.3 

39.9 

37.6 

10* 

148 

91.9 

89.7 

2.5x  10* 

358 

189 

185 

ter  mines  the  minimum  value  of  the  Rayleigh  number  R 
for  which  there  is  steady  finite  amplitude  convection. 
His  result  is 


flmin  =  (vT 3T+  VW, 

(40) 

or 

rm.n=  2pW2  +  T2p, 

(41) 

where  p  nS/R*!3  and  rn'lln  =  (R  m,„  -R*)/r*.  The  value 
of  that  we  compute  on  the  basis  of  Stuart’s  method 
is 

=  2p1/J(l  -  t2)1  n  +  pr2  -  1  .  (42) 


In  the  asymptotic  analysis  of  this  study,  p  —  1  and  t 
are  small.  In  the  limit  p  -  1,  r—  0,  Veronis’  formula 
(41)  predicts  r'mul  =  2,  while  formula  (42)  gives  r'mv>  =  1. 
The  latter  value  is  the  asymptotically  correct  one. 

We  compare  the  predictions  (41)  and  (42)  of  the  semi- 
quantitative  theories  with  the  numerical  work  of  Hup- 
pert  and  Moore.  From  their  solutions,  they  estimate 
R'mtn  as  a  function  of  S  for  fixed  values  of  t.  Table  I 
compares  the  numerical  and  semiquantitative  results 
in  the  case  r  =  10"l/2. 
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